#include "kMST_ILP.h"
#include <sstream>
#include <stack>
#include <set>

using namespace std;

kMST_ILP::kMST_ILP( Instance& _instance, string _model_type, int _k ) :
	instance( _instance ), model_type( _model_type ), k( _k )
{
	n = instance.n_nodes;
	m = instance.n_edges;
	if( k == 0 ) k = n;
}

void kMST_ILP::solve()
{
	try {
		// initialize CPLEX
		env = IloEnv();
		model = IloModel( env );

		// add model-specific constraints
		if( model_type == "scf" ) modelSCF();
		else if( model_type == "mcf" ) modelMCF2();
		else if( model_type == "mtz" ) modelMTZ();
		else {
			cerr << "No existing model chosen\n";
			exit( -1 );
		}

		// build model
		cplex = IloCplex( model );
		// export model to a text file
		cplex.exportModel( "model.lp" );
		// set parameters
		setCPLEXParameters();

		// solve model
		cout << "Calling CPLEX solve ...\n";
		cplex.solve();
		cout << "CPLEX finished.\n\n";
		cout << "CPLEX status: " << cplex.getStatus() << "\n";
		cout << "Branch-and-Bound nodes: " << cplex.getNnodes() << "\n";
		cout << "Objective value: " << cplex.getObjValue() << "\n";
		cout << "CPU time: " << Tools::CPUtime() << "\n\n";

		bool feasible = isSolutionFeasible();
		if( model_type == "scf" ) cout << "Mode = SCF "<< endl;
		else if( model_type == "mcf" ) cout << "Mode = MCF "<< endl;
		else if( model_type == "mtz" ) cout << "Mode = MTZ "<< endl;
		cout << "Found feasible solution: " << (feasible ? "yes" : "no") << endl;
	}
	catch( IloException& e ) {
		cerr << "kMST_ILP: exception " << e << "\n";
		exit( -1 );
	}
	catch( ... ) {
		cerr << "kMST_ILP: unknown exception.\n";
		exit( -1 );
	}
}

bool isOne(IloNum num)
{
	return 0.999 < num && num < 1.001;
}

bool isZero(IloNum num)
{
	return -0.001 < num && num < 0.001;
}

bool equal(IloNum num, unsigned num2)
{
	return abs(num-(double)num2) < 0.001;
}

bool kMST_ILP::isSolutionFeasible()
{
	if (model_type == "mtz" || model_type == "scf" || model_type == "mcf")
	{
		if (model_type == "mtz")
		{
			for (unsigned i=1;i<n;i++)
			{
				IloNum value = cplex.getValue(*orderVars[i]);
				//cout << "o[" << i << "] = " << value << endl;
			}

			/*for (unsigned j=1; j<n; j++) // j is end vertex of edge
			{
			for (unsigned id : instance.incidentEdges[j])
			{
			Instance::Edge& edge = instance.edges[id];
			unsigned i = (edge.v1 == j) ? edge.v2 : edge.v1; // i is start vertex of edge

			IloBoolVar* edgeVar = adj[n*i + j];

			IloNum value = cplex.getValue(*edgeVar);
			cout << "e[" << i << "," << j << "] = " << value << endl;
			}
			}*/
		}
		unsigned n = instance.n_nodes;

		// count edges outgoing from 0 to find root of the tree
		unsigned numEdgesFrom0 = 0;
		unsigned root = -1;
		for (unsigned i=1; i<n; i++)
		{
			IloNum value = cplex.getValue(*adj[i]);
			if (isOne(value))
			{
				numEdgesFrom0++;
				root = i;
			}
		}

		if (numEdgesFrom0 != 1)
		{
			cout << "More than one outgoing edge from node 0!" << endl;
			return false;
		}

		stack<unsigned> nodes;
		nodes.push(root);

		set<unsigned> discoveredNodes;
		discoveredNodes.insert(root);

		unsigned numVerticesReachable = 1;
		unsigned totalCosts = 0;
		//DFS for connectivity check
		while(!nodes.empty())
		{
			unsigned current = nodes.top();
			nodes.pop();

			for (unsigned id : instance.incidentEdges[current])
			{
				Instance::Edge& edge = instance.edges[id];
				unsigned j = (edge.v1 == current) ? edge.v2 : edge.v1;
				if (j != 0)
				{
					IloNum value = cplex.getValue(*adj[current*n + j]);
					if (isOne(value))
					{
						totalCosts += edge.weight;
						nodes.push(j);
						if (discoveredNodes.find(j) != discoveredNodes.end())
						{
							cout << "Discovered cycle in solution!" << endl;
							return false;
						}
						discoveredNodes.insert(j);
						numVerticesReachable++;
					}
				}
			}
		}
		if(numVerticesReachable!=k){
			cout << "Reachable Vertices:" <<numVerticesReachable<<"  but k="<<k<< endl;
		}

		if(!equal(cplex.getObjValue(), totalCosts)){
			cout << "Value :" <<cplex.getObjValue()<<"  but totalCosts="<<totalCosts<< endl;
		}

		return numVerticesReachable == k && equal(cplex.getObjValue(), totalCosts);
	}
	cout << "Reached end of feasibility check" << endl;
	return false;
}

// ----- private methods -----------------------------------------------

void kMST_ILP::setCPLEXParameters()
{
	// print every x-th line of node-log and give more details
	cplex.setParam( IloCplex::MIPInterval,1 );
	//cplex.setParam( IloCplex::MIPDisplay, 2 );
	cplex.setParam( IloCplex::MIPDisplay, 0 );
	//cplex.setParam(IloCplex::RelaxPreInd,0);
	//cplex.setParam(IloCplex::PreInd,0);
	//cplex.setParam(IloCplex::HeurFreq ,-1);
	//cplex.setParam(IloCplex::VarSel ,4);
	//cplex.setParam(IloCplex::PerInd,1);
	// only use a single thread
	cplex.setParam( IloCplex::Threads, 1 );
	cplex.setParam(IloCplex::WorkMem ,10000.0);
}

void kMST_ILP::modelSCF()
{
	// ++++++++++++++++++++++++++++++++++++++++++
	// TODO build single commodity flow model
	// ++++++++++++++++++++++++++++++++++++++++++

	unsigned n = instance.n_nodes;
	unsigned m = instance.n_edges;

	numVariables = 0;
	numConstraints = 0;

	adj = new IloBoolVar*[n*n];
	IloNumVar** adjFlow = new IloNumVar*[n*n];
	vector<IloBoolVar*> allEdgeVars;
	ostringstream oss;


	// add edge variables and flow variables,
	// for each edge between vertices 1,...,n add two directed edge variables
	// for 0, add only outgoing edges
	for (Instance::Edge& e : instance.edges)
	{
		oss << "edges_" << e.v1 << "-" << e.v2;
		IloBoolVar* var1 = new IloBoolVar(env, oss.str().c_str());
		numVariables++;
		adj[e.v1*n + e.v2] = var1;
		allEdgeVars.push_back(var1);

		if (e.v1 != 0)
		{
			// we don't need outgoing flow from vertex 0 as variable
			IloNumVar* flow1 = new IloNumVar(env);
			adjFlow[e.v1*n + e.v2] = flow1;

			oss = ostringstream();
			oss << "edges_" << e.v2 << "-" << e.v1;
			IloBoolVar* var2 = new IloBoolVar(env, oss.str().c_str());
			numVariables++;
			adj[e.v2*n + e.v1] = var2;
			allEdgeVars.push_back(var2);

			IloNumVar* flow2 = new IloNumVar(env);
			adjFlow[e.v2*n + e.v1] = flow2;
		}
	}

	// only one edge outgoing from 0
	IloExpr firstEdge(env);
	for (unsigned i=1; i<n; i++)
	{
		firstEdge += *adj[i];
	}
	model.add(firstEdge == 1);
	numConstraints++;

	// only flow on selected edges, flow must not be larger than k-1 (except for flow from 0)
	// no flow on unselected edges
	for (unsigned i=1; i<n; i++)
	{
		for (unsigned id : instance.incidentEdges[i])
		{
			Instance::Edge& edge = instance.edges[id];
			unsigned neighbor = (edge.v1 == i) ? edge.v2 : edge.v1;
			if (neighbor != 0)
			{
				IloNumVar* flowVar = adjFlow[i*n+neighbor];
				IloBoolVar* edgeVar = adj[i*n+neighbor];

				// no flow on unselected edges, at most k-1 flow on selected edges
				model.add((k-1) * *edgeVar >= *flowVar);

				// at least 1 flow on selected edges
				model.add(*flowVar >= *edgeVar);
				numConstraints+=2;
			}
		}
	}

	// at most one incoming edge per vertex
	for (unsigned i=1; i<n; i++)
	{
		IloExpr incomingEdges(env);
		for (unsigned id : instance.incidentEdges[i])
		{
			Instance::Edge& edge = instance.edges[id];
			unsigned neighbor = (edge.v1 == i) ? edge.v2 : edge.v1;
			incomingEdges += *adj[neighbor*n + i];
		}
		// at most one incoming edge per vertex
		model.add(incomingEdges <= 1);
		incomingEdges.end();
		numConstraints++;
	}

	// incoming flow - outgoing flow = incoming edges
	for (unsigned i=1; i<n; i++)
	{
		IloExpr incomingFlow(env);
		IloExpr outgoingFlow(env);
		IloExpr incomingEdges(env);

		// vertex 0 sends out k flow
		incomingFlow += *adj[i] * k;

		for(unsigned id : instance.incidentEdges[i])
		{
			Instance::Edge& edge = instance.edges[id];
			unsigned neighbor = (edge.v1 == i) ? edge.v2 : edge.v1;
			incomingEdges += *adj[neighbor*n+i];
			if (neighbor != 0)
			{
				incomingFlow += *adjFlow[neighbor*n+i];
				outgoingFlow += *adjFlow[i*n+neighbor];
			}
		}

		model.add(incomingFlow - outgoingFlow == incomingEdges);
		numConstraints++;

		incomingFlow.end();
		outgoingFlow.end();
		incomingEdges.end();
	}

	// objective
	IloExpr obj(env);
	for (Instance::Edge& e : instance.edges)
	{
		obj += *adj[e.v1*n + e.v2] * e.weight;
		if (e.v1 != 0)
		{
			obj += *adj[e.v2*n + e.v1] * e.weight;
		}
	}
	model.add(IloMinimize(env, obj));
	numConstraints++;

	delete[] adjFlow;
	cout << "Num Variables: " << numVariables << endl;
	cout << "Num Constraints: " << numConstraints << endl;
}

EdgeVar* createEdgeVar(IloEnv& env, unsigned start, unsigned end, unsigned weight, unsigned n)
{
	EdgeVar* edgeVar = new EdgeVar();

	ostringstream oss;
	oss << "edge_" << start << "->" << end;
	edgeVar->e = new IloBoolVar(env, oss.str().c_str());

	edgeVar->start = start;
	edgeVar->end = end;
	edgeVar->weight = weight;

	edgeVar->f.resize(n);

	for (unsigned i=1;i<n;i++)
	{
		oss = ostringstream();
		oss << "f_" << start << "->" << end << "_(" << i << ")";
		edgeVar->f[i] = new IloNumVar(env, oss.str().c_str());
	}

	return edgeVar;
}

void kMST_ILP::modelMCF()
{
	// ++++++++++++++++++++++++++++++++++++++++++
	// TODO build multi commodity flow model
	// ++++++++++++++++++++++++++++++++++++++++++

	unsigned n = instance.n_nodes;
	unsigned m = instance.n_edges;

	adj = new IloBoolVar*[n*n];

	incoming.resize(n);
	outgoing.resize(n);

	// create EdgeVar instances, and fill directed incidence lists
	for (Instance::Edge& e : instance.edges)
	{
		EdgeVar* edgeVar = createEdgeVar(env, e.v1, e.v2, e.weight, n);
		incoming[e.v2].push_back(edgeVar);
		outgoing[e.v1].push_back(edgeVar);
		adj[e.v1 * n + e.v2] = edgeVar->e;
		numVariables++;

		if (e.v1 != 0)
		{
			edgeVar = createEdgeVar(env, e.v2, e.v1, e.weight, n);
			incoming[e.v1].push_back(edgeVar);
			outgoing[e.v2].push_back(edgeVar);
			adj[e.v2 * n + e.v1] = edgeVar->e;
			numVariables++;
		}
	}

	// v0 sends out k commodities in total
	IloExpr sourceOutflow(env);
	for (EdgeVar* edgeVar : outgoing[0])
	{
		for (unsigned c=1; c<n; c++)
		{
			sourceOutflow += *(edgeVar->f[c]);
		}
	}
	model.add(k == sourceOutflow);
	sourceOutflow.end();
	numConstraints++;

	// v0 has 1 outgoing edge
	IloExpr expr(env);
	for(EdgeVar* edgeVar : outgoing[0])
	{
		expr += *(edgeVar->e);
	}
	model.add(1 == expr);
	expr.end();
	numConstraints++;

	// each flow is between 0 and the value of the corresponding edge variable
	// selected edges have at least one flow
	for (auto& list : outgoing)
	{
		for (EdgeVar* edgeVar : list)
		{
			for (unsigned c=1; c<n; c++)
			{
				model.add(0 <= *(edgeVar->f[c]));
				model.add(*(edgeVar->f[c]) <= *(edgeVar->e));
				numConstraints+=2;
			}
		}
	}

	// for all nodes: number of incoming edges equals in-flow of its commodity.
	// I.e. if there is an incomging edge on node k, exactly 1 of commodity k most flow in.
	// If there is no inflow on node k, there is no in-flow of commodity k.
	for (unsigned c=1; c<n; c++)
	{
		IloExpr sumOfFlows(env);
		IloExpr sumOfIncomingEdges(env);

		for (EdgeVar* edgeVar : incoming[c])
		{
			sumOfFlows += *(edgeVar->f[c]);
			sumOfIncomingEdges += *(edgeVar->e);
		}

		model.add(sumOfFlows == sumOfIncomingEdges);
		numConstraints++;
		sumOfFlows.end();
		sumOfIncomingEdges.end();
	}

	// on all nodes, for all commodities except its own: preserve commodity
	for(unsigned i=1; i<n; i++)
	{
		for (unsigned c=1; c<n; c++)
		{
			if (i == c)
				continue;

			IloExpr inflow(env);
			IloExpr outflow(env);

			for(EdgeVar* edgeVar : incoming[i])
			{
				inflow += *(edgeVar->f[c]);
			}
			for (EdgeVar* edgeVar : outgoing[i])
			{
				outflow += *(edgeVar->f[c]);
			}

			model.add(inflow == outflow);
			numConstraints++;
		}
	}

	// each node has at most one incoming edge
	for(unsigned i=1; i<n; i++)
	{
		IloExpr incomingEdges(env);
		for (EdgeVar* edgeVar : incoming[i])
		{
			incomingEdges += *(edgeVar->e);
		}
		model.add(incomingEdges <= 1);
		numConstraints++;
	}

	// minimize SUM(weight_i*x_i)
	modelObjectiveFunction();
}

//Kevin Streichers Version of MCF
//Kevin Streichers Version of MCF
void kMST_ILP::modelMCF2()
{
	// ++++++++++++++++++++++++++++++++++++++++++
	// TODO build multi commodity flow model
	// ++++++++++++++++++++++++++++++++++++++++++

	//Number of Vertices
	unsigned n = instance.n_nodes;
	//Number of Edges
	unsigned m = instance.n_edges;

	numVariables = 0;
	numConstraints = 0;

	//2D Adjacency Matrix adj = (a_{ij}) storing pointers to IloBoolVar if there exists an edge from i to j and nullptr else
	adj = new IloBoolVar*[n*n];

	incoming.resize(n);
	outgoing.resize(n);

	// create EdgeVar instances, and fill directed incidence lists
	// EdgeVar is a a struct storing edges in the format:
	// Edge e_{ij} goes from i to j with weight w and there is a flow_k for each commodity k
	// f[0] is nullptr, for 1 <= i <= n-1: f[i] is flow of commodity i on this edge
	for (Instance::Edge& e : instance.edges)
	{
		//Each undirected edge but edges from artificial source node 0 are split into 2 arcs
		//Arc Forward
		EdgeVar* edgeVar = createEdgeVar(env, e.v1, e.v2, e.weight, n);		
		incoming[e.v2].push_back(edgeVar);
		outgoing[e.v1].push_back(edgeVar);
		adj[e.v1 * n + e.v2] = edgeVar->e;
		numVariables++;
		//Arc Back
		if (e.v1 != 0)
		{
			edgeVar = createEdgeVar(env, e.v2, e.v1, e.weight, n);
			incoming[e.v1].push_back(edgeVar);
			outgoing[e.v2].push_back(edgeVar);
			adj[e.v2 * n + e.v1] = edgeVar->e;
			numVariables++;
		}
	}
	//CONSTRAINT 1 - The source sends out k packets
	IloExpr flowOutOfSource(env);
	for (EdgeVar* edgeVar : outgoing[0])
	{	
		for (unsigned c=1; c<n; c++)
		{
			flowOutOfSource += (*edgeVar->f[c]);
		}
	}
	//RELAXTION <= k seems to give a little more speed.
	model.add(flowOutOfSource <= k);
	numConstraints++;
	flowOutOfSource.end();
	cout << "Num Constraints (1): " << numConstraints << endl;

	//CONSTRAINT 2 - Source has 1 outgoing edge
	// - - explanation
	// - - Removing the source has no effect on the kmst anymore
	IloExpr edgesOfSourceSelected(env);
	for(EdgeVar* edge:outgoing[0]){
		edgesOfSourceSelected += *edge->e;
	}
	model.add(edgesOfSourceSelected == 1);
	numConstraints++;
	edgesOfSourceSelected.end();
	cout << "Num Constraints (2): " << numConstraints << endl;

	//CONSTRAINT 3 - Number of Edges is |E| = |V|-1 = (k-1) + 1 from artificial source
	IloExpr numberOfEdges(env);
	for (auto& list:outgoing)
	{		
		for(EdgeVar* edge:list){
			numberOfEdges += *edge->e;
		}
	}
	model.add(numberOfEdges == k);
	numConstraints++;
	numberOfEdges.end();
	cout << "Num Constraints (3): " << numConstraints << endl;

	//CONSTRAINT 4 - f_ij <= |V|-1*x_e = k*x_e
	//- - RAIDL 3.8
	//- - The flow on an edge can not exceed k or 0 if the edge is not selected. Source edge has k*1 as upper bound for the flow. Each node then consumes 1 (its own commodity)
	//- - RAIDL 3.10
	//- - x <= f[c] < k*x
	//- - if	  x=0 : 0 <= f <= 0  -> f=0
	//- - else if x=1 : 1 <= f <= k -> 0 <= f <=k

	//- - it turned out be be faster to use many small constraints for each of the flows this is why i do not use k*x but just 0<=f[c]<=x for each commodity
	for (auto& list:outgoing)
	{			
		for(EdgeVar* edge:list){
			IloExpr flow(env);
			for(unsigned c=1;c<n;c++){
				flow += *edge->f[c];
				model.add(*edge->f[c] <= *edge->e);
				model.add(0 <= *edge->f[c]);
			}			
			//model.add(flow <= k * (*edge->e));
			//model.add((*edge->e) <= flow);
			numConstraints += 2;
			flow.end();
		}
	}
	cout << "Num Constraints (4): " << numConstraints << endl;


	//CONSTRAINT 5
	//- - 1) only 1 incoming arc per node -> tree if connected!
	//- - 2) MCF constraints
	for(unsigned i=1; i<n; i++)
	{//For each node
		IloExpr selectedIncomingEdges(env);
		for(EdgeVar* edgeVar : incoming[i])
		{
			selectedIncomingEdges += *(edgeVar->e);
		}
		numConstraints++;
		//1) See 5.1)
		model.add(selectedIncomingEdges <= 1);


		//2) See 5.1)
		for (unsigned c=1; c<n; c++)
		{//check if all commdodities are passed by except the own.
			IloExpr inflow(env);
			IloExpr outflow(env);
			for(EdgeVar* edgeVar : incoming[i])
			{
				inflow += *(edgeVar->f[c]);
			}
			for (EdgeVar* edgeVar : outgoing[i])
			{
				outflow += *(edgeVar->f[c]);					
			}

			if (i == c){				
				// 0 <= selectedIncomingEdges <= 1
				// 0 -> we have not selected this node, it is already ensured that we can not have a flow here, but we need this
				// to adapt the special case of a k-mst where we can have nodes which do not recieve a packet because they are not selected.
				// 1 -> we have selected this node and therefore we also need to consume one of the packets
				model.add(inflow - outflow == selectedIncomingEdges);
			}else{	
				model.add(inflow - outflow == 0);
			}
			numConstraints++;
			outflow.end();
			inflow.end();
		}
		selectedIncomingEdges.end();
	}

	cout << "Num Constraints (5): " << numConstraints << endl;

	modelObjectiveFunction();

	cout << "Num Variables: " << numVariables << endl;
	cout << "Num Constraints: " << numConstraints << endl;
}

void kMST_ILP::modelMTZ()
{
	// ++++++++++++++++++++++++++++++++++++++++++
	// TODO build Miller-Tucker-Zemlin model
	// ++++++++++++++++++++++++++++++++++++++++++

	unsigned n = instance.n_nodes;
	unsigned m = instance.n_edges;

	unsigned numVariables = 0;
	numConstraints = 0;

	adj = new IloBoolVar*[n*n];
	vector<IloBoolVar*> allEdgeVars;
	ostringstream oss;


	// add edge variables,
	// for each edge between vertices 1,...,n add two directed edge variables
	// for 0, add only outgoing edges
	for (Instance::Edge& e : instance.edges)
	{
		oss << "edges_" << e.v1 << "-" << e.v2;
		IloBoolVar* var1 = new IloBoolVar(env, oss.str().c_str());
		numVariables++;
		adj[e.v1*n + e.v2] = var1;
		allEdgeVars.push_back(var1);

		if (e.v1 != 0)
		{
			oss = ostringstream();
			oss << "edges_" << e.v2 << "-" << e.v1;
			IloBoolVar* var2 = new IloBoolVar(env, oss.str().c_str());
			numVariables++;
			adj[e.v2*n + e.v1] = var2;
			allEdgeVars.push_back(var2);
		}
	}

	// add vertex order variables and constrain to 1 <= var <= k
	orderVars.push_back(nullptr);

	for (unsigned i=1; i<n; i++)
	{
		oss = ostringstream();
		oss << "order_" << i;
		IloNumVar* var = new IloNumVar(env, oss.str().c_str());
		numVariables++;
		orderVars.push_back(var);

		model.add(1 <= *var);
		model.add(*var <= k+1);

		numConstraints+=2;
	}

	// only one edge outgoing from 0
	IloExpr firstEdge(env);
	for (unsigned i=1; i<n; i++)
	{
		firstEdge += *adj[i];
	}
	model.add(firstEdge == 1);
	firstEdge.end();

	numConstraints++;

	// at most one incoming edge per vertex
	// a vertex can only have outgoing edges if it has an incoming edge
	for (unsigned i=1; i<n; i++)
	{
		IloExpr outgoingEdges(env);
		IloExpr incomingEdges(env);
		for (unsigned id : instance.incidentEdges[i])
		{
			Instance::Edge& edge = instance.edges[id];
			unsigned neighbor = (edge.v1 == i) ? edge.v2 : edge.v1;
			incomingEdges += *adj[neighbor*n + i];
			if (neighbor != 0)
			{
				outgoingEdges += *adj[i*n + neighbor];
			}
		}
		model.add(incomingEdges <= 1);
		model.add(incomingEdges * (k-1) >= outgoingEdges);

		// if there are no incomingEdges, order is k+1
		model.add(*orderVars[i] >= (k+1) * (1-incomingEdges));

		outgoingEdges.end();
		incomingEdges.end();

		numConstraints += 2;
	}

	// ordering of vertices in tree (distance to root), u_0 := 0
	// for edges outgoing from vertex 0
	for (unsigned i=1; i<n; i++)
	{
		model.add(*adj[i] <= *orderVars[i]);

		numConstraints++;
	}

	// all other edges
	for (unsigned j=1; j<n; j++) // j is end vertex of edge
	{
		for (unsigned id : instance.incidentEdges[j])
		{
			Instance::Edge& edge = instance.edges[id];
			unsigned i = (edge.v1 == j) ? edge.v2 : edge.v1; // i is start vertex of edge
			if (i != 0) // we already considered edges outgoing from 0
			{
				IloBoolVar* edgeVar = adj[n*i + j];
				IloNumVar* orderStart = orderVars[i];
				IloNumVar* orderEnd = orderVars[j];

				model.add(*orderStart + *edgeVar <= *orderEnd + ceil(k/2.0) * (1-*edgeVar));
				numConstraints++;
			}
		}
	}

	// select k edges in total (including the first edge)
	IloExpr totalEdges(env);
	for (IloBoolVar* var : allEdgeVars)
	{
		totalEdges += *var;
	}
	model.add(totalEdges == k);
	numConstraints++;

	// objective
	IloExpr obj(env);
	for (Instance::Edge& e : instance.edges)
	{
		obj += *adj[e.v1*n + e.v2] * e.weight;
		if (e.v1 != 0)
		{
			obj += *adj[e.v2*n + e.v1] * e.weight;
		}
	}
	model.add(IloMinimize(env, obj));
	numConstraints++;

	cout << "Num Variables: " << numVariables << endl;
	cout << "Num Constraints: " << numConstraints << endl;

}

//TODO: Check if this function is applicable for all 3 methods
void kMST_ILP::modelObjectiveFunction()
{
	IloExpr costs(env);
	for (auto& list : outgoing)
	{
		for (EdgeVar* edgeVar : list)
		{
			costs += *(edgeVar->e) * (double)edgeVar->weight;
		}
	}
	model.add(IloMinimize(env, costs));
	numConstraints++;
	costs.end();
}


kMST_ILP::~kMST_ILP()
{
	// free global CPLEX resources
	cplex.end();
	model.end();
	env.end();
}
